A Mean-field Calculation for the 
Three-Dimensional Holstein Model 1 

Chuan Liu a and Qiang Luo a ' b 

a School of Physics 
Peking University 
Beijing, 100871, P. R. China 

b State Key Laboratory for Mesoscopic Physics 
Peking University 
Beijing, 100811, P. R. China 

Abstract 

A path integral representation appropriate for further Monte Carlo simulations is 
derived for the electron-phonon Holstein model in three spatial dimensions. The 
model is studied within mean-field theory. Charge density wave and superconducting 
phase transitions are discussed. 
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1 Introduction 



The BaPh\- x Bi x O^ alloys [1] are notable ancestors to the high temperature 
superconducting cuprates [2,3]. These compounds exhibit superconductivity 
in the Pb-rich composition range 0.05 < x < 0.3, with a maximum T c w 13K 
near x ~ 0.25. A metal-semiconductor transition is observed at x ~ 0.35, 
with the semiconducting behavior continuing to the end-member compound 
BaBi0 3 [1]. 

By substitutional K doping at the A(Ba) site, Ba x K\_ x BiOz [4] has been 
found to exhibit superconductivity with T c ps 30K. Band-structure calcula- 
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tions [4,5] indicate that the conduction bands in both materials consist of 
o"-antibonding combinations of (Pb,K)-Bi(6s) and 0(2p) states. These bands 
are gradually filled in rigid-band fashion with increasing x until they are half 
filled in BaBi0 3 . Here, a combination of Fermi-surface nesting and the strong 
coupling of the conduction-band states near E F to bond-stretching O displace- 
ments leads to a commensurate charge-density-wave (CDW) order, thereby 
opening a semiconducting gap at Ep, explaining the semiconducting behav- 
ior of BaBi0 3 . The conduction bands may be well represented by simple 
tight-banding models. In contrast to the cuprate high temperature super- 
conductors, these materials are three-dimensional and copper-free, rendering 
magnetic mechanisms inapplicable. Furthermore, conduction electrons couple 
strongly to breathing type oxygen phonons over a wide composition range. 

It appears to us from this picture that the three-dimensional Holstein Hamil- 
tonian on a cubic lattice might constitute a microscopic model for these ma- 
terials. The Holstein Hamiltonian is: 



<xy>,(r x x 

Here, c^ a creates an electron of spin a at site x, and b^ is the creation operator 
for a local phonon mode of bare frequency lu located at site x. The electrons 
have a one-electron overlap t between near-neighbor sites on a cubic lattice, 
ji is the chemical potential and g is the electron-phonon coupling constant. 

For g = 0, the free electrons in the three-dimensional Holstein model have a 
tight binding band structure: e(k) = — 2t(cos k\ + cos/c2 + cos k 3 ). It has per- 
fect nesting only at half-filling. As noted above, this band structure describes 
BaPbi- x Bi x 03 and Ba x K\- x Bi0 3 quite well. The Einstein phonon in the 
Holstein model corresponds to the oxygen breathing mode phonon in real ma- 
terials. When x increases,the rigid-band fashion filling of BaPb\^. x Bi x 3 and 
Ba x K\_ x Bi0 3 corresponds to the variation of the chemical potential fi in the 
Holstein model which determines band filling. All these features indicate the 
connections between the Holstein model and these real materials. It should be 
noted that the boson in the Holstein model could also represent an "exciton", 
a "plasmon", or any other intermediate bosons. The Holstein model therefore 
has a wider range of applicabilities. 

The idea that a single Hamiltonian (1) describes both the semiconducting 
and superconducting phases of BaPb\_ x Bi x Oz and Ba x K\_ x BiOz is very ap- 
pealing. The basic electron-electron interaction arising from the exchange of 
phonons is given by 



V ^ = ~ ,2 -772 ( 2 ) 
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with u m = 2mirT. In the limit wo/t ^> 1 with g 2 /u>o finite, the Holstein 
model goes over to a negative-U Hubbard model with an effective U equal 
to —2g 2 /u [6]. Whereas the basic electron-electron interaction in a Hubbard 
model is instantaneous and local, as is characterized by a constant U, that in 
the Holstein model is retarded, i.e., V(iu m ) depends on frequency (see Eq. (2)). 
This retardation effect will play an important role in subsequent discussions. 

While reports on investigations for the three-dimensional Holstein model have 
been rare, the two-dimensional Holstein model has been previously investi- 
gated by a number of techniques, including quantum Monte Carlo simulations 
[6]. Some of the conclusions from these studies are expected to hold in three 
dimensions. According to these studies, two competing types of correlations 
exist in the Holstein model, s-wave superconducting and Peierls CDW. At 
and near half-filling, CDW long range order prevails, further away from half- 
filling, s-wave superconductivity dominates. As u /t increases, the boundary 
separating CDW and superconducting phases moves toward half-filling. In the 
negative U Hubbard model, which is a limiting case of the Holstein model, 
CDW can only be stabilized at half-filling. It is interesting to note that as 
u /t is reduced and retardation becomes important, a region of filling around 
the half-filled band is dominated by CDW correlations. It is reasonable to 
assume that a CDW phase is insulating. We thus expect that, unlike the Hub- 
bard model, an insulating (semiconducting) phase can be stabilized away from 
half-filling due to retardation in the Holstein model. 

As pointed out by Mattheiss and Hamann [5], the mechanism by which the 
semiconducting properties of BaPh\- x Bi x Oz and Ba x K\^ x BiOj, are extended 
over intermediate range of compositions is less well understood. The semi- 
conducting behavior at half-filling, i.e., BaBiO^, may be understood from a 
combination of Fermi-surface nesting and strong coupling to breathing mode 
oxygen phonons. Both ingredients are necessary to open a semiconductor 
gap at the Fermi surface. The gap should disappear quickly with the loss 
of Fermi surface nesting away from hall-filling. Weber [7] proposed that this 
semiconducting behavior was due to a combination of static incommensurate 
"breathing-type" charge density waves and chemical "ordering waves" on the 
Pb-Bi sublattice. This mechanism, however, has some difficulties with band- 
structure calculations and perturbative arguments [5]. 

It thus appears from these discussions that retardation in the Holstein model 
might be responsible for the stabilization of the semiconducting behavior away 
from half-filling. Ba x Ki^ x BiO^ simply has a phonon frequency higher than 
that of BaPbi^ x Bi x 3 . Consequently, the metallic regime of these alloys is 
extended closer to half-filling, where the deformation potentials for breathing- 
type O displacements and the electron-phonon interaction are at maximum, 
thus resulting in a higher T c . As superconductivity occurs near a CDW tran- 
sition, the coherence length is expected to be short, in agreement with exper- 
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imental observations. 



Several previous studies [8,9] have used certain types of negative-?/ Hubbard 
models to describe BaPbi_ x Bi x 3 systems. The CDWs in Hubbard and Hoi- 
st ein models differ in one respect. CDW in the former is complete in the sense 
that alternating lattice sites are doubly occupied and empty; while in the 
latter, charge is only partly transferred with an amount set by retardation. 
Infrared [11] and x-ray-photoemission-spectroscopy [12] measurements suggest 
only marginal differences in the charge distributions at the Bi sites [5] , lending 
further support to the Holstein model description of these materials. 

Symmetry arguments [10] can be used to classify electronic phase transitions 
in the three-dimensional Holstein model. At half-filling, superconducting and 
CDW channels are degenerate. In this case the two-component superconduct- 
ing order parameter (complex scalar) and the one-component CDW order pa- 
rameter (real scalar) may be considered jointly as a three-component vector 
order parameter. Away from half-filling the degeneracy between the supercon- 
ducting and CDW phases is broken. Thus phase transitions of CDW at half- 
filling, away from half-filling, and superconductivity belong to the universality 
classes of the three-dimensional Heisenberg, Ising, and XY models respectively. 
The corresponding phase transitions in BaPb\_ x Bi x 3 and Ba x K\_ x Bi0 3 are 
thus identified. It will be interesting to explore consequences of this classifica- 
tion in these systems. 



2 Path integral representation of the model 



To obtain a path integral representation for the (grand) canonical partition 
function S = TV exp[— /3H], one first expresses the bosonic creation and an- 
nihilation operators in terms of a bosonic field $ x and the corresponding 
conjugate momenta n x defined as: 




It is easy to see that $ x and n x defined above satisfy the usual commutation 
relations: [ $ x , n y ] = i<5 xy . The Hamiltonian of the model now becomes: 



H = -tJ2 c L( c x+A<r + c*-iur) - ^[fi , v / 2^$x + lA n 



+ Eif + y^. (4) 
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The inverse temperature /3 = 1/T is then divided into N t slices in the "tem- 
poral" direction, each of size a = f3/N t . Inserting the completeness relation 
for the fermionic coherent states and the bosonic variables at any given time 
slice labeled by r: 



xLcr J\xcr / 

X I $r, X ; VV,x CT )(^r+l,x CT ; $t,x |= 1 • (5) 



the partition function can be transformed into: 

N t -1 

n 

T=0 



, N t -1 

„ — n 



where the bosonic and fermionic integrals are over all field components: <3> TX , 
Vv,x<r an d Vv.xo- with r = 0, 1, • • • , N t — 1. In the temporal direction, periodic 
and anti-periodic boundary conditions are understood for the field $ and 
respectively. 

The matrix element in the above expression can be evaluated easily using the 
properties of the coherent states and the fact that the temporal lattice spacing 
a ~ 0. It is convenient to introduce a re-scaled field A x via: 

A x = g>/2u $ x , (6) 



where we have used the notation x to label the four- dimensional site (r, x), 
and g is the dimensionless coupling measured by t, i.e. g = g/t. Similarly, we 
will use the notation: fi = ji/t and u)o = uo/t. Finally, the partition function 
of the model is given by [13,14]: 



S= / VA exp f~E^ ( ^° 2q2) *% ) det(^[A] 2 ) , 

3 

M[A] xy = atY(5 x+ly + 5 x _ ly ) + [1 + at(fi + A,)]^ - C x A-b, y , (7) 
i=\ 

where we have used x = (x, s) to collectively label 3 + 1 dimensional lattice 
sites. The sign function ( Xo — ± is present to insure the anti-periodic boundary 
condition for the fermion fields in the temporal direction. 

It is noted that the real fermion matrix M. is not Hermitean. However, since the 
matrix Ai is real and the up and down spin electrons share the same fermion 
matrix, the effective action is positive definite for any value of the chemical 
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potential. Therefore, this path integral representation for the Holstein model 
is suitable for Monte Carlo simulations [13,14] without the infamous "sign 
problem" . 



3 The charge density wave phase 



To perform a mean-field analysis for the possible charge density wave long 
range order, we assume that the effective action: 

Seff[A x ] = J2A ^~ d [ + 2 Ujy)xy A - 2Tr\ogM[A x \ , (8) 

x,y 4 5 



has a minimum at the "ferrimagnetic" background field A X MF ^ = M + (— ) X A, 
characterized by the magnetization parameter M and the staggered magneti- 
zation parameter A. In this background field, the fermion matrix Ai, which 
will be denoted as Ai^°\ has a simple block diagonal form in four- momentum 
space: 

M pq [A (MF) ] = M$ = [1 - e-*° + at(jl - e p )]5 pq + atA5 q , p+Q , (9) 



where p and q are two four-momenta, the temporal components of which being 
the usual fermionic Matsubara frequencies. The parameters jl and e p designate 
fi + M and e p /t, respectively, and Q equals (0, n, n, ir). As we will show shortly, 
a non-vanishing staggered magnetization order parameter A signals a CDW 
phase in the system while the value of M represents the "renormalization 
effect" of the chemical potential. 

We decompose the scalar field A x as: 

A x = Ai MF ^ +A X = M+ (-) X A + Ax , 



and expand the effective action (8) around the mean- field solution Ai MF \ The 
fact that A^ x MpS) being a saddle point of the effective action insures that the 
term linear in A x is absent in this expansion. Therefore, to second order in 
A x , we obtain: 



S eff [A] = S^[A (MF) ] + si%[A] , 
S (0) [Ai M F)] = ^oa(M 2 + A 2 ) _ M {MF) 

sW[A x ]=Tr(MW- 1 6MM {0) - 1 5M) . (10) 
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This serves as our starting point for the mean-field calculation. In this paper, 
only the lowest order mean-field results will be studied. These results are 
contained in S^L. 

Since the mean-field fermion matrix Ai^ only couples momentum p with 
momentum p + Q, the lowest order effective action = — logS^) can be 
easily obtained: 



S ( %[A {MF) ] = QuJ ° a ( M2 ± A ' 



J^log (l + atjl - e- ipo ) 2 -aH 2 E v 



where we have introduced the energy spectrum in the CDW phase: 

E 2 p = ei + A' . (12) 



We therefore see that, if a non-zero solution for the order parameter A exists, 
the energy spectrum develops an energy gap characterized by A. 

It turns out that the summation over the temporal component p® in eq. (11) 
can be performed analytically. Therefore, we obtain: 



S. 



(0) _ VLuo Q a(M 2 + A 2 ) 

- ^ log 1 + (l + at(j2 - 4 



Nt 



1 + (l + at(fl + E p ] 



Mi 



• (13) 



In most cases, one is only concerned with the so-called continuum limit: a — ► 
with N t a = f3 being fixed. Then, 



(l + at(fl + E p fj * ~ exp (pt{Jl + £ F 



The free energy per lattice site: / = — logH/(/5V3), apart from irrelevant 
constant terms, is thus obtained Q to the lowest order in mean-field theory in 
the continuum limit: 

/(°) = ^(M 2 + A 2 ) - -L ^ log[ooah + cosh(^E p )] . (14) 



2 We use V3 to designate the number of sites of the three dimensional lattice. 
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The chemical potential fi is related to the total electron number via: (N) = 
<91ogS/ ({3d fj,). Using the doping fraction x = (N) /V 3 — 1 and taking the 
continuum limit, this equation reads: 

x - 1 y sinh(ffl/x) 

V 3 V cosh(fStjl) + cosh(f3tE p ) 



It is seen that the solution fi will have the same sign as x, in particular, at 
half-filling (x = 0), the solution to Eq. (15) is: ft = 0. 

The requirement that the background field A^ 1 ^ being a minimum of the 
effective action can be utilized to derive the so-called gap equations. The 
equation for M yields the solution: 

M — 2(1 + x)g 2 /uo , (16) 



where 1 + x is given by eq. (15) , while the equation for A 2 in the continuum 
limit reads: 

<^o_ = J_ y smh((3tE p )/E p 

W V 3 V cosh(/ft/2) + cosh(ptE p ) ' 



At a given doping level x, one needs to use both Eq. (15) and Eq. (17) to 
solve for ft, and A. The critical temperature T c = l/f3 c for the CDW phase 
is obtained from Eq. (17) and Eq. (15) by setting the order parameter A to 
zero, or equivalently, by setting E p = e p for a given value of x: 

J_ sinh(/j c te p )/e p 
2^ 2 ~ U 3 V cosh((3 c tft) + cosh(/5 c te p ) ' 1 ' 



It is also instructive to calculate the ensemble average of the local electron 
number operator n x = J2 a c lo- c xo-- To lowest order in mean- field theory, it is 
given by: 

(0)-l 



(nx> = 2 + 2(X^-; +1)X ) . 



Keeping a finite and after some algebra, one obtains: 

W = ^ + + o(«) = ^ [M + H-A] + 0(«) , (19) 



As we claimed before, a non- vanishing order parameter A opens up an energy 
gap in the spectrum and signals a charge density wave order in the system. 
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Another useful quantity to characterize the CDW phase transition is the CDW 
susceptibility Xcdw defined by: 



Pxcdw = f E (-r y K,xn , y > . (20) 

3 x,y,r 



By virtue of translational invariance, one obtains: 



m X ,y 



Since the fermion propagator for up and down spins are identical, we have: 



Xcdw = E eW*^(x)$(uMy)$(x)) . (21) 



This expression can be represented in momentum space by a loop diagram. 
We will calculate the susceptibility in the phase where A = for simplicity. 
To lowest order in mean-field theory, the susceptibility is given by: 

X<b) ° w = ^ [1 + at(fi - lp) - e-™] [1 + at(fi + e p ) - F^] ' (22) 



where we have used the fact that: e p = — e p+ Q. Using similar techniques as 
before, the summation over po can be performed analytically, and in the con- 
tinuum limit, the above equation is reduced to: 

(o) J_ v / J_\ sinh(/?e p ) m] 

Xcdw ^ ^ £ J cosh(/%/i) + cosh (^ p ) " ^ ^ 



A better approximation, the Random Phase Approximation (RPA), for the 
CDW susceptibility can be obtained by summing over a geometric series of 
X^c with one phonon exchange. In this approximation, the CDW susceptibility 
becomes: 

(o) 

(RPA) Xcdw (0A \ 

Xcdw ~ (n n.~ 2 i * \ (o) ' ^ > 

1 - (2fitg 2 u )xcdw 



which becomes divergent exactly at the critical temperature T c , as is evident 
by comparing with Eq. (18). 



9 



4 The superconducting phase transition 



In the previous section, we have discussed the possibility of the CDW phase 
transition in the Holstein model. This phase transition is characterized by a 
long range order parameter A which is diagonal in electron spins. It turns 
out that the system also permits a phase with long range order which is off- 
diagonal in the space of the electron spins. This phase is the conventional 
BCS-type superconducting phase with Cooper pairing. To see such a phase 
transition, the phonon field is integrated out from the path integral represen- 
tation of the partition function. Then, another set of boson fields, namely a 
real scalar field $ x , a complex scalar field A^ and its complex conjugate A*, 
are integrated in. Introducing the four component fermion field notation: 



tlx 
tlx 
\ *P2x J 



(25) 



the action of the model can be written as: 



S = -y ax V[A,<S>] ax ^p y 



Ag 2 u a 



% ■ 



(26) 



xy 



where the index a and f3 run from 1 to 4 and the anti-symmetric fermion 
matrix V[A, ^} aX ;f3 y is given by: 



'0 -atA*J xy 



-A4°f - 




















(27) 



where the matrix Ai^ is given by: 

M<$ = atj:(5 x+ly + 5 x _ ly ) + (1 + atjl)5 xy - ■ (28) 

i=i 



The effective action of the system is obtained by integrating out the fermion 
fields, resulting in the Pfaffian (the square root of the determinant) of the 
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corresponding fermion matrix T>[A, $]: 

= f E I A.I' H- (=|±^) ^ * s - lTr W , „ . , 29 ) 

We now proceed to search for a constant mean-field solution of the type: 
A x = Abcs and $ x = $ which minimizes the effective action. In this case, the 
fermion matrix V is diagonal in four-momentum space, i.e. V pq = V^\p)5 pq 
and it is given by: 



P(°)(p) = 






a ~t^*BCS 


-M<®(p)* 













-M(0)( p )* 










atA BC s 







atA BCS 






with .M(°)(p) = 1 + a£(/t — e p ) — e tpo . Therefore, to lowest order in mean-field 
theory, the effective action is found to be: 



-I^log[|l-e-^ + at(/i-e p )| 2 + (a0 2 |A BC5 | 2 ] . (31) 

Since the fermion contribution is independent of $o, it is immediately seen that 
the solution for <3> is $o — 0- The summation over Matsubara frequencies can 
be performed analytically. In the continuum limit the result is: 

where the spectrum E p is given by: 



E p = ^fi-e p y + ^A BCS . (33) 

Chemical potential \i is determined from band filling. Using the doping fraction 
parameter x, the chemical potential satisfies the following equation: 

It is noted that, at half-filling x = and the solution is /i — due to the 
symmetry of e p . The superconducting energy gap A BC s is obtained from the 



11 



famous BCS gap equation: 



1 X ^ 



9 2 V 3 ^\E r 



(35) 



The critical temperature for the superconducting phase transition is obtained 
by setting the BCS gap parameter A BC s to zero in the above equation: 



g 2 



l 

v 3 



E 



t 



|/i-e r 



tanh ( 



(36) 



The fermion propagator is given by the inverse matrix of T> which again is 
diagonal in momentum space. If we denote: 



" p 



P> (0) U(p) , 



(37) 



then the matrix [D^] x (p) is given by: 





( ° 


atA BC s 


M {0) {j>Y 





1 


-atA BCS 










D(p) 










—atA* BCS 




\ o 


-M^ip) 


a ^BCS 






(38) 



with D{p) = \M (0 \p)\ 2 + \atA BCS \ 2 . 

We can also calculate the superconducting susceptibility defined by: 



Xbcs 



x,y 



(39) 



For simplicity, we will evaluate the superconducting susceptibility in the phase 
with vanishing A B cs (high temperature phase). After some algebra, we obtain 
to the lowest order: 



(o) - \ - 
Xbcs ~ 77 2^ 



1 



1 



tanh ( ^7 £pl 



^3 p P\H ~ e P 
Under RPA, the susceptibility becomes: 



(40) 



ARPA) 
Xbcs 



(o) 
Xbcs 



1 - iPtgyu ) X Pcs 



(41) 
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which becomes divergent exactly at f3 c determined by Eq. (36). 



5 Results and discussions 




Fig. 1. In the lower half of the figure, the critical temperatures for the CDW 
phase transition are plotted as a function of g 2 /ojq for four different doping lev- 
els: x = 0.00, —0.15, —0.25, —0.35. In the upper half of the figure, the order pa- 
rameters for the CDW phase transition are plotted as a function of the tem- 
perature T/t for a given value of ojq/9 2 = 0.2 for four different doping levels: 
x = 0.0,-0.15,-0.25,-0.35. 

Mean-field equations discussed in previous sections can be evaluated numer- 
ically and the results for the critical temperatures and the order parameters 
are summarized in Fig. 1 and Fig. 2. In the lower half of the figures, we plot 
the critical temperatures as a function of g 2 /uj for four different doping levels: 
x = 0.00, —0.15, —0.25, —0.35. The critical temperatures for the superconduct- 
ing phase transition are represented by the four lines that increase with g 2 /ooo 
near the bottom of Fig. 2. They show rather weak dependence on doping. The 
corresponding lines for the CDW phase transition illustrate a much stronger 
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Fig. 2. In the lower half of the figure, the critical temperatures for the BCS su- 
perconducting phase transition are plotted as a function of g 2 /ujq for four different 
doping levels: x = 0.00,-0.15,-0.25,-0.35. In the upper half of the figure, the 
order parameters for the BCS superconducting phase transition are plotted as a 
function of the temperature T/t for a given value of ujo/g 2 = 0.2 for four different 
doping levels: x = 0.0, -0.15, -0.25, -0.35. 

dependence on doping compared to the superconducting critical temperatures. 
It is noted that, for a given value of g 2 /ooo, the critical temperatures for the 
CDW phase transition decrease as the doping is increased. At exactly half 
filling, the critical temperature for the CDW phase transition is always higher 
than that of the superconducting phase transition. While the issue of the sta- 
bility of the CDW and superconducting phases can only be addressed by a 
more sophisticated calculation, here we simply take the magnitude of the crit- 
ical temperatures and order parameters as a measure of the tendency that the 
system is likely to be in one of these phases. Thus, we expect the system to 
be in a CDW phase at half filling. Away from half filling, the critical tempera- 
ture for the CDW phase transition eventually becomes smaller than that of the 
superconducting phase transition. This is a signal that the CDW phase will 
eventually becomes unstable when the system is doped away from half filling. 
The superconducting phase will become dominant at low temperatures. 
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In the upper half of the figures, we plot the order parameters A = A/t for the 
CDW and the superconducting phase transitions as a function of temperature 
T/t at Coo/g 2 = 0.2 for the same four different doping levels. It is seen that, for a 
given value of u>o/g 2 , both order parameters decrease as the doping is increased. 
However, the order parameter for the CDW phase decreases much more rapidly 
than the superconducting order parameter when the system moves away from 
half filling. These features indicate that the CDW phase is most favored at 
half filling, in agreement with our qualitative expectations discussed in the 
introduction. 



6 Conclusions 

In this paper, we present our mean-field analysis of the three dimensional Hoi- 
st ein model. A path integral representation of the model is obtained which is 
suitable for further Monte Carlo simulations. The critical temperatures and 
the order parameters are obtained within mean-field approximation for both 
the CDW phase and BCS superconducting phase transitions. Our results are 
consistent with the picture that the CDW phase is most favored at half-filling. 
Away from half filling, the usual BCS superconducting phase becomes impor- 
tant. Higher order calculations are necessary to further explore the competition 
between different electronic phase transitions. 
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